rm(list=ls())
library(foreign)
library(gdata)

setwd("/scratch/cds2083")
# Before starting this, combine all of the estimates
# from the different imputation runs (stored in the
# int-out-imp-1.csv through int-out-imp-10.csv spreadsheets 
# into one XLSX spreadsheet called “int-out-imp-all.xls.” 
# An example of what this new spreadsheet should look like
# is provided in the replication archive.
# That file also includes the regression, matching, and 
# IPW estimates produced in Stata.

intres.sc <- read.xls("int-out-imp-all.xls", sheet=1)
intres <- intres.sc[nrow(intres.sc):1,]

allBounds <- c(	intres$bavg+intres$moeavg, 
				intres$bavg-intres$moeavg,
				intres$lbreg,
				intres$ubreg,
				intres$lbmatch,
				intres$ubmatch,
				intres$lbnipw,
				intres$ubnipw)

plotmax <- ceiling(max(allBounds)*10)/10
plotmin <- floor(min(allBounds)*10)/10
pdf(file="int-out-rev.pdf", width=8, height=8)
par(mar=c(5,0,0,0))
plot(c(-1.2, plotmax),
		c(0,nrow(intres)), 
		type="n",
		axes=F,
		xlab=bquote(Delta~Recid.~Index~Popn.~Mean),
		ylab="")
abline(v=0)
points(intres$bavg, 1:nrow(intres), pch=19)
segments(	intres$bavg+intres$moeavg, 
			1:nrow(intres), 
			intres$bavg-intres$moeavg,
			1:nrow(intres))

points(intres$breg, 
			1:nrow(intres)-.15, 
			pch=21)
segments(	intres$lbreg, 
			1:nrow(intres)-.15, 
			intres$ubreg,
			1:nrow(intres)-.15)

points(intres$bmatch, 
			1:nrow(intres)-.3, 
			pch=21)
segments(	intres$lbmatch, 
			1:nrow(intres)-.3, 
			intres$ubmatch,
			1:nrow(intres)-.3)

points(intres$bnipw, 
			1:nrow(intres)-.45, 
			pch=21)
segments(	intres$lbnipw, 
			1:nrow(intres)-.45, 
			intres$ubnipw,
			1:nrow(intres)-.45)

text(rep(-1.2, nrow(intres)), 1:nrow(intres), 
			intres$intname,
			pos=4)

text(rep(plotmin, nrow(intres)), 1:nrow(intres), "Ensemble IPW", cex=.75)
text(rep(plotmin, nrow(intres)), (1:nrow(intres))-.15, "WLS", cex=.75)
text(rep(plotmin, nrow(intres)), (1:nrow(intres))-.30, "Matching", cex=.75)
text(rep(plotmin, nrow(intres)), (1:nrow(intres))-.45, "Naive IPW", cex=.75)

axis(1, round(10*(seq(plotmin,
			plotmax,
			by=.1)))/10,
			cex=.75)
dev.off()